home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1995 April / macformat-023.iso / Shareware City / Developers / RLaB / toolbox / expm.r < prev    next >
Encoding:
Text File  |  1994-02-21  |  1.1 KB  |  66 lines  |  [TEXT/RLAB]
  1. //-------------------------------------------------------------------//
  2.  
  3. //  Syntax:    expm ( A )
  4.  
  5. //  Description:
  6.  
  7. //  The expm function computes the matrix exponential, or exp (A),
  8. //  where A is a square matrix.
  9.  
  10. //  Reference: "Matrix Computations" Gene  Golub, and Charles Van
  11. //          Loan. Section 11.3.
  12.  
  13. //-------------------------------------------------------------------//
  14.  
  15. expm = function ( a_ )
  16. {
  17.   local (X, a, c, cX, d, k, N, n, q, j);
  18.  
  19.   if (a_.nr != a_.nc) { error ("expm: A must be square"); }
  20.  
  21.   //
  22.   // Scale A if neccesary, such that the infinity-norm of A divided
  23.   // by some power of 2 is less than 1/2.
  24.   //
  25.  
  26.   j = norm(a_, "i");
  27.   if (j > 0) 
  28.   { 
  29.     j = max ([0, fix (log (j)/log (2))+2]);
  30.   }
  31.   a = a_/2^j;
  32.   n = a.nr;
  33.  
  34.   //
  35.   // k = 1
  36.   //
  37.  
  38.   c = 1/2;
  39.   X = a; 
  40.   d = eye (n, n) - c*a;
  41.   N = eye (n, n) + c*a;
  42.   q = 6;
  43.  
  44.   for (k in 2:q)
  45.   {
  46.     c = c*(q - k + 1)/(k*(2*q - k + 1));
  47.     X = a*X;
  48.     cX = c*X;
  49.     N = N + cX;
  50.     d = d + (-1)^k*cX;
  51.   }
  52.  
  53.   //
  54.   // Solve d*F = N for F, overwrite N for efficiency
  55.   //
  56.  
  57.   N = d\N;
  58.  
  59.   for (k in 1:j)
  60.   {
  61.     N = N*N;
  62.   }
  63.  
  64.   return N;
  65. };
  66.